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Progress Report 

The major accomplishment during this contract period has been the incor- 
poration of the multiblock derivative interface algorithm into the flux-split 
Euler code developed here at MSU. The code can now be used to compute 
flow on any multiblock configuration provided the block boundary surfaces 
coincide. There is no need for any type of grid line or slope continuity at 
the block boundaries. We do not permit the overlapping of blocks so that 
all interpolation can be done on surfaces, which is a two- dimensional in- 
terpolation problem. Since our Euler code is a cell-centered scheme, there 
arose some unexpected problems as to how the phantom points should be up- 
dated. This is mainly a problem with the second-order scheme where there 
are two surrounding layers of phantom points that must be updated after 
each time step. After trying several alternatives, it was decided that linear 
extrapolation from the values on the block boundary worked best. These 
problems did not arise with the MacCormack method that was used in our 
earlier developmental work. One of our concerns in the earlier work with 
this interface procedure was whether the updating at the boundary could 
be done efficiently. The update procedure is an implicit algorithm and the 
boundary values are computed at each time step using an ADI scheme. In 
a three-dimensional test case, it was observed that the number of ADI it- 
erations has practically no effect on the numerical solution. Therefore, we 
are using only one ADI iteration per time step. Thus, in terms of operation 
count, the update procedure is competitive with any explicit procedure. It 
should be noted that these conclusions are preliminary since we only recently 
began computations on three- dimensional multiblock grids. 

The following example is included to illustrate the type of results that 
have been seen. The geometry is a simple three- dimensional axisymmetric 
body with a spherical grid system. A supersonic flow is computed with a zero 
angle of attack. Discontinuities in grid lines occurs when the grid is refined 
near the leading and trailing edges. The contours appearing in the following 
figure were plotted from a solution computed using the derivative interface 
procedure. As a check on the procedure, a solution was also computed by 
extending the grid lines into the neighboring block and interpolating solution 
values. There was no noticeable difference in the two solutions. 

The presentation ” Derivative Interface Conditions for Multiblock Grids” 
made at the Third International Conference on Numerical Grid Generation in 
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Figure 1: Interface procedured applied to a three-dimensional grid with grid 
line discontinuity 

Barcelona included our results as of June 1991. The paper which appeared in 
the Proceedings was submitted earlier in the spring and was attached to our 
last semiannual report. We have included with this report the paper entitled 
” Linear Variational Methods and Adaptive Grids”. The paper was recently 
accepted for publication in Computers and Mathematics with Applications . 
The manuscript was submitted over a year ago and contains work that was 
done about two years ago. However, since there is now a renewed interest in 
variational methods, the work is still timely. 



To appear in Computers and Mathematics with Applications 


Linear Variational Methods and Adaptive 
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1 Introduction 

The progression of larger and larger supercomputers has not eliminated 
the need for adaptive grids in computational fluid dynamics. Instead, it 
has brought demands for more detailed resolution in the flow field and the 
solution of more complicated fluid flow problems. Regardless of the chosen 
algorithm, the resolution of shock waves and boundary layers demands ex- 
tremely close spacing of grid points. Shock waves can be especially difficult 
to deal with because they appear in the interior of the computational region 
and may be a transient phenomena. 

The adaptive grids considered in this report will be structured grids 
and will be constructed by moving existing grid points into regions where 
solution variables have large derivatives. In recent years the term adaptive 
grid has expanded to encompass such topics as local grid refinements and 
fine grid overlays. These techniques can be used without having to deal with 
problems of grid distortion, but they do require an additional level of data 
structure and cannot be easily incorporated into an existing grid generation 
/ flow solver package. 

The redistribution of grid points in the construction of adaptive grids 
can be done using either direct algebraic methods or iterative methods that 
are derived from variational problems or elliptic boundary value problems. 
Algebraic methods are the simplest and generally involve a spline fit of the 
points based on some equidistribution principle. As with all algebraic meth- 
ods, there is no way of insuring against grid folding (negative Jacobians), 
especially if the redistribution is done in more than one coordinate direc- 
tion. Iterative methods for constructing adaptive grids have now been used 
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routinely for several years. The initial work on the variational based meth- 
ods was done by Brackbill and Saltzman [3]. The foundations for the most 
popular elliptic method were laid by Anderson and Steinbrenner [1] using 
the system of elliptic partial differential equations developed by Thompson 
et al [9]. Presently the equations of Anderson and Steinbrenner are most 
widely used for three-dimensional problems since the Euler equations for the 
variational method are very lengthy. 

Comprehensive review articles on adaptive grids have been written by 
Eiseman [5] and Thompson [8]. The elliptic grid generation equations of 
Thompson are accepted as the best system for constructing arbitrary curvi- 
linear grids. There is only one problem that may arise in implementing the 
method. The elliptic equations which must be solved numerically are nonlin- 
ear and iterative methods for solving the equations do not always converge. 
This is particularly true in the case of adaptive grid construction where 
there are large variations in grid spacing. The nonlinearity also implies that 
numerical algorithms will only be locally convergent and, therefore, con- 
vergence will only occur if the iterative sequence starts with a reasonable 
estimate of the final grid. While none of these problems are insurmountable, 
it has lead to an interest in grid generation methods based on the solution 
of linear systems of algebraic equations. Linear algebraic equations may 
be derived from the discretization of linear systems of partial differential 
equations or from the direct discretization of certain variational problems. 
Most recently work in this area has been reported by Castillo et al [4]. A 
somewhat similar approach was followed by Kennon and Dulikravich [6]. 
The latter method was based on a discrete optimization problem while the 
former considered a continuous variational problem. Depending on the dis- 
cretization used in the variational problem, the methods may or may not 
generate the same system of linear algebraic equations. A comparison of 
grids generated by linear and nonlinear methods can be found in the reports 
of Castillo et al and Mastin et al [7]. 

The objective of this report is to derive linear systems for grid gener- 
ation which naturally tend to generate nonfolding grids. The key to the 
method is to attempt to preserve the aspect ratios of the grid cells when 
transforming between the physical and computational regions. Surprisingly, 
this fact was noted many years ago by Barfield [2]. Only two-dimensional 
grids were considered, in which case the appropriate generating equations 
can be derived from properties of conformal mappings. The same equa- 
tions will be derived here in a slightly different way which leads to grid 
generating equations in three dimensions and also equations for construct- 


2 



ing two and three-dimensional adaptive grids. Examples have demonstrated 
that nonfolding grids can be constructed for many nonconvex regions where 
methods that use smoothing formulas derived from Laplace’s equation fail. 
The adaptive grids constructed using linear methods are similar to those 
constructed by nonlinear methods and are more orthogonal in some cases. 
However, the methods which are presented here should not be considered as 
a replacement for existing elliptic and variational schemes. There were cases 
where they failed to generate a satisfactory grid while the existing nonlinear 
schemes were successful. On the other hand, they are still attractive because 
all the traditional iterative schemes to solve the linear systems of algebraic 
equations will always converge regardless of the initial grid coordinates. 

2 Two-Dimensional Grid Generation 

The variational approach will be used to generate the difference equa- 
tions. This leads directly to the generalization to three dimensions and the 
equations for adaptive grids. Since one of the problems in the construction 
of adaptive grids is skewness, the development will begin with the variational 
principle of conformal mappings. 

The adaptive grid will be constructed from the mappings in Figure 1. 
An arbitrary bounded simply-connected region R in the xy-plane can be 
conformally mapped onto a rectangular region with four given points of R 
mapping to the vertices of the rectangle. If one of the mapping functions 
is scaled, the resulting transformation maps R onto a square S. The aspect 
ratio of the rectangular region is a conformal invariant of R called its module 
and denoted by M. The transformation from the square onto R minimizes 
the following integral 

JJ s M ( x l + vl) + Jj( x l + yl) d l* dl '- (!) 

The conformal mapping also satisfies boundary orthogonality constraints 
which will not be considered in this development. One-dimensional stretch- 
ing transformations can be composed with conformal mappings without de- 
stroying the one-to-one and orthogonality properties of the mapping. Thus, 
one-dimensional transformations defined by the equations 

M = «(0> 0 < £ < 1 

v — b(rj), 0 < 77 < 1 
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Figure 1: Composite Mappings for a Two-Dimensional Grid 


will be included in the mapping from the square onto the region R. The 
resulting transformation minimizes the integral 



i q'(0 

MV(Tj) 


{*1 + Vi) ^ d V- 


( 2 ) 


One of the objectives of an adaptive grid is to equidistribute some quantity 
over the total region. The transformations a and b can be chosen from an 
equidistribution property. Suppose that /(£) is a positive function defined 
for 0 < £ < 1. A one-dimensional mapping will equidistribute the function 
/ on the unit interval if 

f(O d H = cd£, 


where c is a constant to be determined. The equation can also be written 
as 



( 3 ) 


Integrating from £ = 0 to £ = 1 gives 
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or 


= Uo W) dr ) 


-i 


f( r ) 

Therefore, an equidistribution of /(£) in the fit /- plane can be accomplished 
from the one- dimensional stretching transformation given by the above dif- 
ferential equation or explicitly as 



Since the mapping from the 7^-plane to the zy-plane is generated from 
properties of conformal mappings, the resulting grid will also locally exhibit 
an approximate equidistribution property. An equidistribution in the other 
coordinate direction would be done by choosing the function 6(77) in a similar 
fashion. If orthogonality is essential, then the adaptivity functions are re- 
stricted to one-dimension. However, since the orthogonality condition is not 
really our aim, two-dimensional equidistribution functions will be allowed 
so that the integral in (2) can be generalized. Two functions /(£, 77) and 
g(£, 77) will be considered with / controlling adaptation in the ^-direction 
and g controlling adaptation in the 77-direction. It will be assumed that 
these functions have been normalized along grid lines so that the constant 
in (3) is at least approximately c = 1. The adaptive grid is constructed by 
minimizing the integral 
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(x 2 r, + y 2 r,)dtdTl. 


(4) 


Two things are still needed to make this an easily solvable problem. The 
boundary correspondence must be given. Equidistribution of a weight func- 
tion along a curve is a simple algebraic procedure and it will be assumed 
that this has been done before the interior grid points are to be calculated. 
No attempt will be made to duplicate the boundary correspondence of the 
conformal mapping since that would lead to a nonlinear problem. The quan- 
tity M must also be calculated. Computational experience has shown that 
only a crude approximation is necessary. If this grid generation scheme is 
to be an iterative scheme, the initial grid may be an algebraic grid on R. 
In that case one could compute an average aspect ratio of all cells of the 
algebraic grid. Using the mapping notation, we would compute M as 


M = 




d£dr). 
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where 


* = 

y - y(Z,v) 


is any mapping of the square region onto the region R. A simpler procedure 
which avoids the need for an initial algebraic grid is to define M as the ratio 
of the average lengths of opposite sides of R. 

A discretization of the integral (4) leads to a least-squares problem which 
in turn yields linear equations for the coordinates x and y. Since both 
equations are the same, we let r = ( x,y ), i and j denote the grid point 
indices, and write the system as 


~ ?i,j) + Pi-'jft- i.i " ?ij) 

+Qi,j+ 1 ( f iJ+ 1 - rij) + Q tJ _ i (rij- 1 - fij ) = 0 

where 


P *±|.i 


Qij±h 




M 




Since P and Q are positive, this linear system of equations for the grid 
points ( ) is diagonally dominant and can be solved using any direct 
of iterative method. 

It should be noted that this grid adaptation method causes grid cluster- 
ing by controlling the cell aspect ratio. The same control function cannot 
be used to cluster simultaneously in the £ and r) directions because if / = g, 
the integral (4) reduces to integral (1). Of course, grid adaptation could be 
achieved by introducing a weight function in (1), that is, replace dfidu by 
w(fj.,i , )dfxdu, but then the method would loose much of the orthogonality 
inherent in its development. 


3 Three-Dimensional Grid Generation 

The basic concept in the two-dimensional development is preservation of 
cell aspect ratios in conformal mappings. Although conformal mappings do 
not extend to three dimensions, this basic concept can be formulated. The 
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Figure 2: Composite Mappings for a Three-Dimensional Grid 

initial assumption is that given a three-dimensional region R, a one-to-one 
mapping of a rectangular region Q onto R can be found which minimizes 
the integral 

/// IM 2 + ||r„|| 2 + \\r w \\ 2 dudvdw 

provided the length, L, width, W, and height, H, of Q axe properly defined 
and suitable boundary conditions are imposed. Although this may seem a 
plausible statement, no results of this type have been found in the literature 
and it is doubtful that it is true for all regions. Thus the three-dimensional 
generalization begins on a less firm footing and the value of the final grid 
generation algorithm must be judged by its ability to generate nonfolding 
grids of practical interest. 

A sequence of mappings will be examined as illustrated in Figure 2. The 
first step will be the mapping of the rectangular region Q onto a unit cube, 
C. This can be done by scaling so that 

u = L/jl, v = Wv, w = Hu. 

The integral (5) becomes 

J I Jc + if *'<*■' (5) 
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In view of our basic assumption it is now evident that the mapping of the 
unit cube onto an arbitrary region should depend on the dimensions of the 
region. Here we have three parameters to approximate instead of one as in 
integral ( 1 ). Proceeding as in two dimensions, if an initial grid is given on 
R from say an algebraic mapping of the cube, then the following ratios can 
be computed. 


Mi = 


WH 

~L~ 

LH 


■IIL-i 


dfidudu 


*-S-//jW*** 


Once again the geometric ratios of the region R are approximated by com- 
puting the average ratios of all the grid cells in R. Grid adaptivity can also 
be included by first considering one-dimensional stretching and then gen- 
eralizing the control functions so that the degree of stretching is allowed 
to vary between grid lines. If the stretching transformations map the unit 
cube in £77^- space to the unit cube in g,uu> -space and the functions to be 
approximately equi distributed along the £, 77, and C grid lines are /, g, and 
h, then the integral (6) to be minimized becomes 


/ I Sc Ml ^ l|fd|2 + A j|fj2 + M3 7F l|fd|2 * dT]dc - 

Three-dimensional adaptivity can also be achieved by weighting the differ- 
ential as described in the previous section. A direct discretization of the 
variational problem leads to the following diagonally dominant linear sys- 
tem. 


where 


^«'+ ~ l,j,k ~ Ti'j'k) 

~ ~ f i,j,k ) 


P 

Q 

R 


A77AC/ 
A £gh 
A£ A( g 
A77 fh 
Af A77 h 

“ac 77 


Mi 


M 2 


m 3 
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Conclusions 


The objective of this report has been to present a reliable linear system 
for grid generation in two and three dimensions. The method is robust in 
the sense that convergence is guaranteed but is not as reliable as other non- 
linear elliptic methods in generating nonfolding grids. The construction of 
nonfolding grids depends on having reasonable approximations of cell aspect 
ratios and an appropriate distribution of grid points on the boundary of the 
region. Some guidelines have been included on approximating the aspect 
ratios, but we can offer little help on setting up the boundary grid other 
than to say that in two dimensions the boundary correspondence should be 
close to that generated by a conformal mapping. 

It has been assumed that the functions which control the grid distribu- 
tion depend only on the computational variables, £, 77 , and £, and not on the 
physical variables, x,y, and z. Whether this is actually the case depends 
on how the grid is constructed. In a dynamic adaptive procedure where the 
grid is constructed in the process of solving a fluid flow problem, the grid 
is usually updated at fixed iteration counts using the current value of the 
control function. Since the control function is not being updated during the 
iteration of the grid equations, the grid construction is a linear procedure. 
However, in the case of a static adaptive procedure where a trial solution 
is computed and used to construct an adaptive grid, the control functions 
may be recomputed at every step.of the adaptive grid iteration based on the 
current location of the grid points. This latter case gives rise to a nonlin- 
ear system and convergence cannot be guaranteed regardless of the elliptic 
system used for grid generation. 
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